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ABSTRACT 


The  flow  fields  of  a  steady  and  a  periodically  pulsed  two- 
dimensional  turbulent  free  jet  have  been  studied  by  solving 
the  thin  shear  layer  equations  by  the  Keller  Box  method  in 
transformed  variable  form.  A  constant  eddy- viscosity 
formulation  was  used  to  model  the  Reynolds  shear  stress 
term.  For  the  steady  jet,  calculations  agree  well  with 
documented  experimental  data.  Computed  results  of  the  un¬ 
steady  jet  indicate  that  the  mean  flow  characteristics 
follow  closely  those  of  the  steady  jet  and  compare  well  with 
available  experimental  data.  For  sufficiently  high  frequency 
and  amplitude  of  pulsation  or  at  large  streamwise  distance, 
significant  unsteady  effects  occur  in  the  instantaneous 
quantities . 
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NOMENCLATURE 


Unless  otherwise  stated,  the  symbols  used  in  the  text 
have  the  following  meanings. 

a  constant  defined  in  Eq.  (15a) 

b  constant  defined  in  Eq.  (ISb) 

c  constant  defined  in  Eq .  (11) 

F  %  peak-to-peak  oscillation  of  centre-line  velocity 

f  stream  function  in  (S,t,n)  coordinates  defined  in 

Eq.  (15b) 

h  nozzle  width 

M  mean  momentum  flux  in  streamwise  direction 

P  *  mean  pressure 

p*  instantaneous  pressure 

Q  mean  mass  flow  at  any  streamwise  station 

QE  mean  mass  flow  at  nozzle  exit 

T  period  of  pulsation 

t  non-dimensional  time  =  Uc^*t*/h 

U  non-dimensional  mean  x-component  velocity  =  U*/U 

U  (y)  non-dimensional  mean  velocity  profile  at  the 
nozzle  exit 

u  non-diemnsional  instantaneous  x-component  velocity 

u*/U  * 
ci 

u'  non-dimensional  x-component  velocity  fluctuation  = 

u*’/U*ci 

V  non-dimensional  mean  y-component  velocity  =  V*/U*c^ 

non-dimensional  instantaneous  y-component  velocity 
v*/U*ci 
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v'  non-dimensional  y-component  velocity  fluctuation 

v*'/U*cl 

x  non-dimensional  streamwise  distance  -  x*/h 

y  non-dimensional  transverse  distance  -  y*/h 

y. ^  value  of  y  at  which  u  =  yUc 

GREEK  SYMBOLS 

A  difference  between  two  quantities 

c  amplitude  of  pulsation 

v  non-dimensional  kinematic  viscosity  =  v*/ (U*c^h) 

eddy-viscosity  defined  in  Eq.  (10) 

veff  Effective  eddy  viscosity  defined  in  Eq.  (13) 

stream  function  in  (s,t,y)  coordinates  defined  in 
Eq.  (15b) 

(b  phase  angle 

$  function  defined  in  Eq.  (28) 

n  transformed  variable  defined  in  Eq .  (15b) 

oj  non-dimensional  angular  frequency  of  pulsation  = 

o)*h/Uci* 

S  non-dimensional  distance  identical  to  x 

So  constant  defined  in  Eq .  (15a) 

SUPERSCRIPT 

*  dimensional  quantities 
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SUBSCRIPTS 

c  centerline 

i  initial 

00  jet  edge 

Overbar  time-average  quantities 
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1.0  Introduction 


The  steady  two-dimensional  turbulent  jet  has  long  been 
the  subject  of  many  theoretical  and  experimental  investiga¬ 
tions  because  of  its  important  role  in  many  different  types 
of  engineering  applications  such  as  fluidic  and  combustion 
systems  and  because  of  its  significance  in  providing  funda¬ 
mental  understanding  to  the  physics  of  flow  mechanisms  such 
as  turbulence  and  vortex  structures.  Owing  to  its  simplicity 
in  configuration,  the  steady  two-dimensional  turbulent  jet 
has  been  studied  in  detail  experimentally  and  is  regarded  as 
one  of  the  most  well-documented  flows  in  the  literature 
where  sufficient  data  are  available  for  many  practical 
engineering  purposes.  However,  despite  such  extensive 
investigations,  for  example  Heskestad  (1965),  Gutmark  and 
Wygnanski  (1976)  and  Everitt  and  Robins  (1978)  ,  considerable 
scatter  is  found  to  exist  between  the  results  of  various 
workers  even  in  the  mean  flow  parameters  such  as  centre-line 
velocity  decay  rate  and  jet  spreading  rate.  Such  discrepan¬ 
cies  have  been  casually  attributed  to  different  effects  such 
as  Reynolds  number,  aspect  ratio,  nozzle  geometry,  initial 
conditions,  upstream  turbulence  intensities  and  the 
uncertainties  involved  in  the  hot  wire  results  in  regions 
where  reversed  flow  may  occur.  Nevertheless,  no  unified 
agreement  on  the  effects  of  such  factors  and  other  flow 
mechanisms  on  the  flow  development  has  been  reached.  A 
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comprehensive  review  and  evaluation  of  the  experiemntal 
data  on  steady  turbulent  jets  was  given  by  Harsha  (1971) 
and  Rodi  (1975). 

The  analytical  solution  of  the  fully-developed  steady 
turbulent  jet  was  first  sought  by  Tollmien  (1926)  followed 
by  Goertler  (1942)  and  Schlichting  (1965) .  With  the  advent 
of  computer  technology  and  the  rapid  development  of  numeri¬ 
cal  techniques,  the  st  idy  turbulent  jets  have  commonly  been 
used  for  turbulence  modelling  and  as  standard  test  cases  of 
turbulence  models.  Many  numerical  predictions  of  steady 
turbulent  jets  using  various  turbulence  models  have  been 
attempted  such  as  Rodi  and  Spalding  (1970) ,  Launder  et  al 
(1972)  and  Chen  and  Nikitopoulos  (1979)  to  calculate  the 
flow  properties  and  to  complement  the  experimental  results. 

Although  almost  inevitably  unsteadiness  of  varying 
degrees  occurs  in  practice  either  desirably  to  achieve 
certain  favorable  characteristics  or  undersirably  due  to 
the  fluctuations  in  the  surrounding  fluid,  very  few  results 
appear  to  have  been  reported  on  unsteady  turbulent  jets  both 
theoretically  and  experimentally.  Only  until  recently, 
because  of  the  growing  realization  of  the  fundamental  and 
practical  implications  of  an  improved  understanding  of 
unsteady  effects,  the  excitation  of  turbulent  jets  by 
acoustic  (Fiedler  and  Korschelt  (1979))  mechanical  (Bremhorst 
and  Harcn  (1979))  and  fluidic  (Piatt  and  Viets  (1979))  means 
has  received  considerable  attention.  On  the  other  hcnd, 
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closed  form  solutions  of  the  unsteady  jets  which  adequately 
describe  the  flow  development  can  hardly  be  obtained  despite 
the  efforts  of  Pai  (1965)  and  McCormack  et  al  (1966). 
Numerical  solution  of  unsteady  laminar  jets  was  obtained  by 
Kent  (1973).  Although  turbulence  models  have  been  developed 
to  give  sufficiently  accurate  predictions  of  a  wide  variety 
of  steady  flows,  the  applicability  of  such  turbulence  models 
to  unsteady  flows  is  uncertain. 

The  objectives  of  this  study  are  to  apply  the  trans¬ 
formation  developed  by  Lai  and  Simmons  (1978)  to  compute 
steady  and  unsteady  turbulent  jets,  to  add  to  the  under¬ 
standing  of  the  steady  and  unsteady  two-dimensional 
turbulent  free  jets,  to  investigate  the  validity  of  quasi¬ 
steady  approximations  and  to  evaluate  the  suitability  of 
using  turbulence  models  established  for  steady  flows  in 
unsteady  flow  calculations.  The  unsteady  jet  considered 
consists  of  an  initially  steady,  two-dimensional,  turbulent 
free  jet  with  a  sinusoidal  mass  flow  variation  superimposed 
on  it  at  the  nozzle  exit.  The  steady-state  oscillatory  flow 
characteristics  at  any  location  downstream  of  the  nozzle  are 
obtained  by  solving  the  thin  shear  layer  equations  in  two 
spatial  (t,n)  and  one  time  (t)  transformed  coordinates. 
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2.0  Governing  Equations 


The  Navier-Stokes  equations  for  two-dimensional 
incompressible  flow  in  tensor  notations  are  given  by  (see, 
e.g. ,  Hinze  (1975)) 


* 


+ 


x  3p*  3ZU.* 

p  +  V*  3x>  & 


with  the  continuity  equation  being  given  by 


3u-  * 
w  ~  0 


3xi 


(2) 


By  applying  Reynolds  decomposition 


ui*  =  U^*  +  u.^*'  ,  etc. 

to  equation  (1)  and  taking  time  average  with  a  time  scale 
large  compared  with  that  of  the  turbulent  motions  but  small 
compared  with  the  periodicity  of  the  flow,  the  time-averaged 
Reynolds  equations  can  then  be  given  by 


9U. *  3U. * 

i  +  fi  *  ^ 

3t*“  uj  3x7* 


3U.  * 

.  p*  +  i_  +  

v  7x7*  3x 
J 


3  U.  * 


(3) 


Consider  an  unsteady,  two-dimensional,  constant  property, 
turbulent  free  jet  issuing  into  a  stationary  medium.  The 
instantaneous  configuration  is  shown  schematically  in  Fig.  1. 
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Assuming  the  thin  shear  layer  approximations  as  in  Cebeci 
and  Bradshaw  (1977),  the  governing  equations  can  be  obtained 
from  equations  (2)  and  (3)  to  yield  in  rectangular 
coordinate  system 


dU*  dV*  n 
Tx*  Sy*  '  u 


9U*  +  U*  +  V*3U 


+  =  -JL.(y*.il L*  .  u*'v*') 

3y*  Ty*1-  Sy* 


The  boundary  conditions  on  the  centre-line  and  at  the 
edge  of  the  jet  can  be  expressed  respectively  by 


t*  >  0 


=  V*  =  0 


y*  =  y0 


U  *  a  0 


The  time-varying  boundary  condition  at  the  nozzle  exit 
is  given  by 

t*  >  0  x  *  =  0  U*  =  U  *(y*)(l  +  esinw*t*)  (6a) 

where  UQ*(y*)  is  the  mean  velocity  profile  at  the  nozzle 
exit,  e  is  the  amplitude  of  pulsation  and  to*  is  the 
angular  frequency  of  pulsation. 

For  all  locations  downstream  of  the  nozzle,  the 
initial  conditions  are  given  by 


t*  a  0 


x  *  _>  0  U*  =  Ui*(?c  *,y*) 
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where  U^*(x*,y*)  is  the  steady  state  solution  of  equation 
(4)  with  the  boundary  condition  at  the  nozzle  exit  given 
by  equation  (6a)  with  t*  *  0  . 
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2 . 1  Governing  Equations  in  Non-Dimensional  Form: 

In  order  to  enable  the  solution  obtained  to  be  valid 
for  a  family  of  velocity  profiles  which  have  the  same 
normalized  shape  at  the  nozzle  exit,  the  following  non- 
dimensional  variables  are  used:  - 


ucx.y.t)  .  vcx.y.t)  = 

vc i  U  . * 

ci 


x  = 


X* 

ir 


U 


t  = 


ci 


0)  = 


U 


-3F-W' 


Cl 


V*  — , — r-  U*’ V*’ 

v  =  U™ = 

C1  Cl 

where  Uci*  is  the  centre-line  velocity  at  the  nozzle  exit 
at  t*  =  0  . 

Equation  (4)  can  now  be  written  as 


9U  9V 
9x  9y 


0 


(7a) 


9U 

9t 


u|2 

9x 


V—  =  _JL(V—  -  u' v' ) 

v9y  dy''  9y  v  J 


(7b) 


subject  to  the  following  boundary  and  initial  conditions 


t  >  0 


y  =  0 

y  s  oo 


t  >  o 


x  =  o 

x  >  0 


9U 

w 


V 


0 


u  =  o 


U  =  UQ(y)(l  +  esinajt) 
U  =  U.  (x,y) 


(8a) 

(8b) 


t 


0 
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(9a) 

(9b) 


2.2  Turbulence  Modelling 

The  Navier-Stokes  equations  constitute  a  complete  set 
of  equations  of  motion  which  can  in  theory  be  solved  to 
yield  a  solution  for  any  laminar  or  turbulent  flow  field. 
However,  turbulence  comprises  a  wide  range  of  length 
scales  bounded  from  above  by  the  dimensions  of  the  flow 
field  and  bounded  from  below  by  the  diffusive  action  of 
molecular  viscosity.  The  resolution  scale  of  the  smaller 
eddies  which  are  responsible  for  the  decay  of  turbulence  is 
too  small  that  it  precludes  the  use  of  any  existing 
computer.  Futhermore,  very  often  only  the  time-averaged 
properties  are  of  interest  in  engineering  applications 
even  if  the  flow  is  time  dependent.  Hence,  the  Navier-Stokes 
equations  are  time-averaged  first  before  being  solved.  This 
avoids  not  only  the  difficulty  in  representing  all  the 
characteristic  turbulence  scales  but  also  unnecessary  compu¬ 
tations  of  transients  if  only  the  time-averaged  properties 
are  required.  The  process  of  time-averaging  such  as  the 
Reynolds  time-averaging  describe d  in  section  2.0  causes 
the  loss  of  certain  information  contained  in  the  original 
equations  and  results  in  more  unknowns  than  the  governing 
equations  through  the  introduction  of  statistical  correla¬ 
tions  of  fluctuating  velocities  such  as  the  u'v'  term 
in  equation  (7b)  which  are  known  as  apparent  Reynolds 
stresses  and  are  responsible  for  the  actual  momentum 
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transfer.  Attempts  to  derive  additional  equations  for 
those  Reynolds  stresses  will  only  result  in  additional 
unknowns.  Thus  the  time- averaging  process  presents  a 
closure  problem  which  is  to  reduce  the  number  of  unkonwns 
to  equal  the  number  of  governing  equations.  In  order  to 
achieve  this,  the  additional  unknown  quantities  must  be 
modelled  or  approximated  in  terms  of  known  quantities  through 
a  set  of  equations  which,  when  solved  with  the  mean-flow 
equations,  simulate  the  actual  flow  situation.  This  process 
is  generally  termed  "turbulence  modelling." 

Turbulence  models  can  broadly  be  classified  into 
"first-order"  models  in  which  the  mean  flow  equations  are 
solved  without  additional  partial  differential  equations 
for  the  velocity  fluctuation  terms  and  "higher-order" 
models  in  which  transport  equations  for  higher-order  velocity 
correlations  are  solved  with  the  mean  flow  equations.  Since 
in  the  "first-order"  models  additional  partial  differential 
equations  are  not  solved  for  the  turbulence  quantities  which 
are  expressed  in  terms  of  known  quantities  through  certain 
algebraic  formulation,  they  are  known  as  zero-equation  models 
and  "higher-order"  models  which  involve  at  least  one  partial 
differential  equation  for  the  turbulence  quantities  can  thus 
be  termed  one-equation  or  two-equation  models,  etc.  A 
description  of  the  various  types  of  turbulence  models  can 
be  found  in  Launder  and  Spalding  (1972)  and  a  comprehensive 
review  of  the  state-of-art  is  given  by  Reynolds  and  Cebeci 
(1978)  . 


Zero-equation  models  which  are  based  on  empirical 
correlations  of  the  extensive  available  experimental  data 
and  mostly  on  the  eddy-viscosity  and  mixing-length  concepts 
such  as  the  Cebeci-Smith  eddy  viscosity  model  (1974)  have 
been  widely  used  and  proved  to  be  very  successful  in  ob¬ 
taining  very  accurate  predictions  of  non-separating  flows. 

For  more  complex  flows  which  involve  recirculation  and 
separation,  the  more  refined  two-equation  k-e  model 
described  by  Launder  and  Spalding  (1976)  is  the  most  well- 
developed  among  all  other  "high-order"  models.  It  has  been 
recognized  that  although  "higher-order"  models  contain  more 
information  and  hence  simulate  the  flow  situation  more 
realistically  than  the  "first-order"  models,  they  are  more 
difficult  to  solve  and  require  generally  an  order-of 
magnitude  more  computing  storage  and  time  while  in  non¬ 
separating  and  simple  flow  situations,  the  first-order 
models  can  yield  predictions  to  the  same  degree  of  accuracy 
as  the  "higher-order"  models.  It  is  in  the  light  of  this 
philosophy  that  in  this  study,  a  constant-eddy  viscosity 
model  due  to  Prandtl  (1942)  is  used. 

By  employing  the  eddy- viscosity  concept  of  Boussinesq, 
the  fluctuating- velocity  correlation  term  u' v'  in 
equation  (7b)  is  related  to  the  mean  velocity  gradient  through 
the  eddy  viscosity  vt  as  follows: 


uTvr 


=  -  v 


3U 

t3y 


(10) 
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The  eddy  viscosity  vt  can  then  be  obtained  from  the 
constant  eddy-viscosity  model  as 

vt(t)  =  cy1/2  Uc(t)  (11) 

where  c  is  a  constant,  and  y^y2  t^ie  jet  ^alf -width. 

In  this  model,  two  constants  and  c 2  are  used  for 

the  free  jet  as  it  emerges  from  developing  to  fully  developed 
form.  Hence  c  is  given  by 

r  c,  if  U_  £  0.95 
c  =  {  1  C  (12) 

t  c 2  if  Uc  <  0.95 


By  defining  an  effective  eddy  viscosity  such 

that 


v 


eff  = 


vt  +  v 


(13) 


and  using  equation  (ID), equation  (7b)  can  then  be  written  as 


3t 


+ 


3  ,  3U.. 

Jy  tveff5yJ 


(14) 


Hence  the  flow  development  of  a  periodically  pulsed 
turbulent  free  jet  is  governed  by  equations  (7a)  and  (14) 
subject  to  the  boundary  and  initial  conditions  of  equations 
(8)  and  (9)  with  the  constant  eddy  viscosity  model  of  (11) 
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2.3  Governing  Equation  in  Transformed  Variable  Form 

It  has  been  shown  by  Lai  and  Simmons  (1978)  that 
tranformations  based  on  the  similarity  solutions  are 
successful  in  reducing  substantially  the  rate  of 
spread  of  an  unsteady  laminar  jet,  hence  enabling 
the  use  of  variable  grid  sizes  over  predetermined 
regions.  In  this  study,  similar  transformations 
are  employed.  A  dimensionless  transverse  dis¬ 
tance  n  and  a  dimensionless  stream  function  ^  are 
defined  by 

n  *  ay/ +  50)  (15a) 

K?,y,t)  -  b(e  +  ?0)1/2f(c,n,t)  (15b) 

where  a,  b  and  ZQ  are  arbitrary  constants  which  can 
be  varied  to  facilitate  computation  and  x  is  renamed  as 
C  . 

The  function  f  automatically  satisfies  the  continuity 
equation  (7a)  and  equation  (14)  can  be  re-written  in  the 
transformed  coordinates  as 


u  *  ?0)'1/2(veff 


£"]’  ♦  (f 'Y  +  ff’ 


2tc  *  i0nr|f 


-  ♦ 

7c 


Cc  ♦  cj1/2  3f, 


(16) 


"aET 


at 


1 


\ 
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where  prime  denotes  differentiation  with  respect  to  n  , 


f'  -  U  +  5  )1/2U/ab  (17) 

o 

and  a  and  b  are  chosen  for  convenience  such  that 
a/b  =  \ 

The  boundary  conditions  in  equation  (8)  become 


t  >  0 


n 

n 


0 


f"  =  0 
f'  =  0 


f*2U*  t0)|f  -  o 


(18a) 

(18b) 


The  initial  conditions  at  5  =  0  in  the  (  n  ,  t)  plane 
is  obtained  by  writing  equation  (9a)  with  equation  (17)  as 

f'  =  f^frOd  +  esincot)  (19a) 

The  initial  conditions  at  t  =  0  in  the  (s,n)  plane 
are  generated  by  the  solutions  U^(5,n)  of  the  following 
steady  jet  equation  obtained  from  equation  (16) 

Cc  +  cQ)  Evef  f  f "] '  +  (f’)2  +  ff"  =  2(c  +  ;0)[f,H' 

-£”||j  (19b) 


The  effective  eddy-viscosity  can  be  obtained  from 
equation  (13)  using  equations  (11)  and  (15)  to  give 

veff(t)  -  cb  (5  ♦  5g)1/2  n1/2fc,(t)  +  v  (20) 

where  c  is  subject  to  the  conditions  specified  in 
equation  (12) . 
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3.0  Method  of  Solution 


Since  the  initial  conditions  in  the  (S,n)  plane  are 
generated  by  solving  the  steady  jet  equation  (19b)  subject 
to  the  boundary  conditions  in  equation  (18) ,  the  flow 
development  of  the  steady  jet  is  first  studied.  For 
computer  programme  development  purpose,  two  separate 
computer  programs  were  written  for  the  steady  and  unsteady 
jet  respectively  with  the  solutions  of  the  steady  jet 
stored  on  disk  for  use  in  the  unsteady  jet  program  .  How¬ 
ever,  if  necessary,  the  steady  jet  program  can  easily  be 
incorporated  into  the  unsteady  one.  Both  equations  (16)  and 
(19b)  are  parabolic  and  can  be  solved  by  a  marching  procedure. 
The  finite  difference  scheme  employed  is  the  Box  Method 
developed  by  Keller  (1970)  ,  and  described  in  detail  by  Cebeci 
and  Bradshaw  (1977).  The  scheme  has  been  applied  successfully 
to  a  wide  range  of  boundary -layer  type  of  flows  including  both 
time-dependent  and  separating  flows  by  Cebeci  and  his  co¬ 
workers  (e.g.,  see  Bradshaw  et  al  (1980)). 
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Equation  (16)  is  rewritten  as  a  system  of  first  order 


parital  differential  equations 
£'  ■  g 


s'  -  q 

u  +  sor1/2  ^ 


eff 


q]  '  +  g2  +  fq  =  2U  +  ?0)  [g||  - 


Cf 

qs<; 


r  +  r  1/2 

*  (1  ,to)  it] 

a5  St J 


(21a) 

(21b) 

(21c) 


Consider  the  net  cube  shown  in  Fig.  2  and  denote  the  grid 
points  by 


0  ?i+l  =  ?i  +  ri  1  =  1>2...,  I 


t  =  0  t_  ®  t_  .  • 

o  n  n-1  n 


+  K  n  =  1,2,.  .  ,  N 


nn  -  0  * 


j  ’  n;-l  *  hj  J  *  T'2'-'- 


"j  '  n. 


where  i,n,j  are  just  sequence  numbers  and  the  variable  net 
spacings  r^,  kn  and  hj  are  completely  arbitrary. 

In  the  study  of  jets,  a  coarse  grid  spacing  hj  can  be  used 
in  the  vicinity  of  the  centre  line  and  the  edge  of  the  jet 
as  the  transverse  (n)  gradients  of  f  are  small  in 
those  regions.  The  quantities  (f,g,q)  at  points  (?^  ,  tn, 
n j )  are  approximated  by  the  grid  functions  (fj’n  , 

gj’n  ,  qj,n).  Hence  by  using  central  differencing  and 
averaging  about  the  mid-point  (?i >  tn,  _ ^)  ,  Eqs .  (21a) 

and  (21b)  can  be  written  in  the  following  finite  difference 


(22a) 


Ef 


(fi,n  .  £i.n)/h 


3  g 


(g-’n  -  -1,n 


1-h 

q1-? 

4J  -h 


(22b) 


where,  for  example,  the  shorthand  notation  g-V  has  been 
used  for  >s(gj’n  +  gj:J) 

Shorthand  notations  are  introduced,  for  example, 

iz  =  Js(f^,n  +  f-’11'1  +  f^’3?  +  f-’^'1)  =  J*(fj'?  +  f^»?-1) 

A  J  J  J  *  A  J“A  J  J  “'l 

fj  ■  3^ Cf j  * n  +  fj'1,n  +  fj’11"1  +  fi"l>n“l)  _  >4(fj>n+  fj234) 

fn  =  +  fj‘1,n  +  £jl?  +  =  W)'.\  + 


£] 


234 


£i,n-l  +  ^i-l.n  +  ^i-l,n-l 


Eq.  (21c)  can  be  approximated  using  central  differencing  at 

C;i-V  Vv  by 

(Bq)j  -  (Bq)j.i  +  hjUg2)j.jj  +  (fq)j_j,  -  atgjVgj.jj  * 

gj :2_1  +  g]l\  -  2§i-i)  -  qj.j5Cfj.j5  + 


ai8j  -3s}  =  Tj-h 


fi  ,n-l  ‘r  \  r  234«i 

•  2£i-l>  '  ' 

_  Ti-l,n-l 


(22c) 


where  B  *  (S  +  z^j  ) "% 
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a 


m 


i-h 

1 


m1  =  2(?  +  50) 


3/2, i-h 


4[U  +  50)j/a] 


ab  k 


n 


Ti-1 ,n-l 

Tj-h 


,  234  ,  234  _  ,  ,  ,  2,  234  ,  234 

CBq)  i  -  (Bq) j  +  {-  (g  ) ^ _h  -  (fq)^.^  + 

a^82  8j-h  '  0lqj-*J  +  ale3' 


$ 

8i 


g1’?'1 
Sj  .% 

ri ,n-l 


234  ,  - 

8j-h  ^i-1 

2  fi-l 


8?  =  g 


,  *n~ 1  - 


j-h 


2  «i-: 


8,  - 


g1'*’11 
sj  -*5 


-  2  g 


n- 1 


In  the  above  equation  (22c),  the  superscripts  i,n  have 
been  dropped  for  simplicity. 

It  has  been  pointed  out  by  Keller  (1978)  that  in 
approximating  non-linear  terms  such  as  (fq)j-jj  equati°n 
(22c)  with  the  Box  scheme,  there  exists  several  choices 
which  should  not  have  serious  effect  on  accuracy  or  stability 
as  long  as  the  proper  centering  is  maintained.  However,  in 
jet  calculations,  it  has  been  found  in  this  study  that  the 
approximation  of  terms  such  as  (fq)  in  the  form  of  Cf)j 
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results  in  better  numerical  stability  at  the  jet  edge. 
The  boundary  conditions  in  Eq .  (18)  become 


qp ,U  =  0  (23a) 

i  n  2m,  ^  2m«  •  , 

fo  =  C1  +  Fp  '  ^£o  (23b) 

g},n  =  0  (23c) 


The  initial  conditions  at  ?  *  0  in  the  (n,t)  plane  in 
Eq  (19a)  are  given  by 


(£')“  »  (f'0)j[l  +  esin  co(nAt)  ]  (24) 


where  At  is  the  temporal  grid  size. 

The  initial  conditions  at  t  *  0  in  the  (c.n)  plane 
correspond  to  thQ  steady  jet  solutions  of  equation  (19b). 
The  finite  difference  equations  can  be  obtained  similar  to 
the  above  derivation  for  the  unsteady  jet  except  terms  in 
equation  (22)  at  t  =  tn are  equal  to  those  at  t  =  tn 
and  <*1  =  0  . 
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Time  Varying  Velocity  Profile  at  the  Nozzle  Exit 
Because  of  symmetry,  only  half  of  the  jet  needs  to  be 
computed.  The  time-varying  velocity  profile  considered  at 
the  nozzle  exit  is  the  commonly  assumed  "top-hat"  profile 
which  is  given  by 

1  +  esin  .w  t  for  0  4  y  4  1/2 
U(y,t)  =  I  (25) 

|  0  for  1/2  <  y 

In  terms  of  the  transformed  coordinates,  equation  (25)  can 
be  rewritten  using  equation  (17)  as 

(1  +  esin  wt)  for  0  <  n  £  (nOT)i 
for  (Oi  i  y 

where  (n^) ^  *  0.5  a/?0  (26) 


Since  a  discontinuity  exists  in  f'  at  n  *  (nm)i  and  will 
cause  computational  problems,  the  top-hat  profile  has  to  be 


approximated  by 


f  = 


U  *  50) 


(n)  (1  +  esin  wt) 


where 

^(n) 


for  0  <  n  < 

—  c L 

nN2)(l  +  2nN)  for  na  <  n  <  (n0O)i 
for  (n.lj  £  n 


nM 


(28) 
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and  na  is  some  point  in  the  interval  [0,  (n^)- ]  . 

2 

The  function  (1  -  n^) (1  +  2n^)  was  obtained  by  matching  the 
velocity  and  velocity  gradient  profiles  at  n  ■  (n  ) .  and 

03  X 

n  =  na  through  the  following  boundary  conditions  to  ensure 
continuity: 


nN  '  0 

-  1  =  0 

(29a) 

nN  *  1 

-  0  4£  -  0 

(29b) 

f  and  f" 

can  be 

obtained  by  respectively  integrating  and 

differentiating  equation  (27)  with  respect 

to  n  to  yield 

f  = 

U  +  c0) 

1/2 

4 (n)  (1  ♦  esin  fi>t) 

(30) 

i. 

ab 

and  f"  = 

(?  +  50) 
a5 

1/2 

<5 2  (1  +  esin  wt) 

(31) 

*0(n)  =  ! 

1 

fCnJi 

-  na]nN(l  -  +  0.5r^3) 

for  0  <_  n  £  na 
for  n a  <  ti  <  (n0O)i 

(32) 

0 

for  (nJi  £  o 

and 

1 

for 

0  <  n  <  na 

$2Cn)  - 

j  ‘6nN  (1 

'  nN)/ [ Cn„) ±  *  nal  for 

—  n  <_CT’Ji 

0 

for 

(Oi  £  n 

(S3) 

The  value  of  n_  can  be  determined  to  match  the 

a 

experimental  top-hat  profile  as  close  as  possible  and  in 
this  study  it  is  chosen  to  be  the  point  such  that  the 
velocity  profile  is  uniform  over  97.5%  of  the  nozzle. 
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•n  mimmmm 


Adopting  s 
a  and  r 

o 

such  that 


Lmilar  procedure  given  in  Lai  and  Simmons  (1978), 
are  given  the  values  0.3  and  3. 75  respectively 
(n.4  =0.04  and  na  =  0.039. 
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3.3  Solution  of  the  Finite  Difference  Equation 
With  (f?  *  ,  g?*1,  q^**)  known  from  the  solution  of 
Eq  C19b)  and  (fj’1^1'1^'1)  specified  by  Eqs.  (27),  (30) 
and  (31),  Eq.(22)  for  1  <  j  <  J  and  Eq.  (23)  yield  an 
implicit  non-linear  algebraic  system  of  3J  +  3  equations. 
This  system  is  linearized  by  Newton's  Method  by  introducing 
the  perturbed  quantities  (6f,<5g,Sq)  to  yield  after  consider¬ 
able  algebra 


6fj  “ 


6gj 


(sl)j 


where 


Csl>J 

(s2)j 

Cs3>j 

(s4)j 

(s5)j 


(Tl)j 


j  _ i  -  7^-C<Sgj  +  <sg± _1)  =  (ri)J.  (34a) 

6Sj-l  *  +  <Sqj  _iD  -  (rj)..!  (34b) 

Sqj  +  (S2)j  6qj.1  ♦  (S3) j  <Sfj  ♦  Cs4)j5fj.1 

+  C*sVg3  +  Cs6)  j  6sj -1  =  (r25j  (34c) 

hi 

=  Bj  +  T-  Cfj  +  afj-Jj  *  aBl) 
hi 

=  B..1+  *  af  j  -Jj  +  a8x) 

b-:  2  34 

=  r-  fqj  +  ac*j-ii + 

’  /  (4j-l  * 

-  hjCSj  '  “gj.i,  -  ^  ‘  Z1’ 

’  hj  (sj-l  '  agj-H  '  T~  ~ 


fj-r  f;  *  hj  8i->> 
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j - 1 

fr2Jj 


gi-l  '  Si  *hi  qj-<i 


•  Tj:!-"-1  *  CBq)  j  _!  -  (Bq)j  -  hj  {gjJ.%  (*!),.„ 

-  “I'sCSj.ij  *  Sj  8j-l)  *  B  Sj  _jj  -  fj-SjBj-li 


>  j-  1 

>!  *j.H  -  -  «1  8j-v 


The  boundary  conditions  are  given  by 


6ffl  =  0 


<5gj  *  0 


The  above  derivation  is  similar  to  that  outlined  in  Appendix 
C  in  Lai  8  Simmons  (1978).  The  linear  system  Eq.  (3.10)  is 
then  solved  very  effectively  by  the  Block  Elimination  Method 
discussed  by  Cebeci  and  Bradshaw  (1977). 

The  linearized  form  of  the  steady  jet  equation  can 
similarly  be  obtained  and  solved  with  quantities  at 
t  =  tnl  equal  to  those  at  t  *  tn  and  *  0  . 


3.4  Convergence  Criterion 


As  the  governing  equations  for  both  the  steady  and  unsteady 
jet  are  parabolic,  they  are  solved  by  marching  along  the 
t-direction  in  the  case  of  the  unsteady  jet  and  along  the 
^-direction  in  the  case  of  the  steady  jet.  At  a  given  stream- 
wise  (?)  station,  the  linearized  system  Eq .  (34)  is  solved 
by  iterating  at  each  t-station  until  some  convergence 
criterion  is  satisfied.  Iterations  are  terminated  at  each 
t-station  if 

| £t (i+l)  -  <  ei  ,  at  n  =  ncon  (36) 

where  the  value  of  is  prescribed; 

f ' ^  and  are  the  i  th  and  (i+l)th 

iterates  of  f'  respectively;  and  ncon  is  some  point 
where  the  convergence  criterion  is  applied.  In  practice,  it 
is  adequate  to  set  e1  to  be  10"3  . 
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■ 

3.5  Criterion  for  the  Spreading  of  Jet 
Since  the  mass  flow  varies  sinusoidally  with  time,  and 
because  of  the  initially  steep  velocity  gradient,  the  jet 
width  will  be  changing  with  time  at  a  given  streamwise  (?) 
station.  A  criterion  must  therefore  be  set  to  determine 
at  each  t-station,  noting  that  •  Here 

and  nf*  ^  are  the  jet  width  for  the  i  th  and 
(i-l)th  iterations.  The  edge  of  the  jet  is  defined  by  the 
following  two  conditions  - 


f ' J-l  1  -  c2 

(37a) 

f"i  1  i  *3 

(37b) 

where  the  values  of  e2  and  £3  are  prescribed  and  J 
denotes  the  point  at  the  jet  boundary.  Experience  indicates 
that  it  is  sufficient  to  choose 

e2  »  10'2  and  £3  =  10*1 


If  the  criteria  set  out  in  Eq.  (37)  are  satisfied  then 


,  (i  +  l) 


(i) 


Otherwise  n^  points  have  to  be  added  so 


that  J  *  J  -  ,  +  n„  and  the  values  of 
new  old  g 


(£j1,n  gj1,n»  q-,n 


) 


for  the  new 


rij  points  are  obtained  as  follows  - 


ft 


j *n  *  Cnj  -  nj  gj‘n  «■  £j,n 


g-,n  =  g},n 


(38a) 

(38b) 
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0 


(38c) 


Bj,n  =  Bj ,n 


(38d) 


The  same  procedure  is  also  applied  to 


qj’H'1, 

i-1 ,n-l 


Bi,n-1  fi-l,n 
Ri-1 ,n-l 

’  bj 


i*l »n  i-l,n 
j  *  4j 


fi »n-l 
j 

gl.T!-1 

1 


cr 1  >  n~  1 

gj 

^i-1 ,n-l 
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3.6  Criterion  for  the  Attainment  of  Steady  State  Solution 
As  the  mass  flow  at  the  nozzle  varies  sinusoidally  with 
time,  the  steady  state  solution  at  any  streamwise  (?)  station 
downstream  must  vary  periodically  and  can  be  expressed  in  the 
form 

f'U.t.n)  =  fJU.t.n)  +  Z  f£(?,n)  sin  (*wt  +  *4)  (39) 

1 

In  general,  for  a  periodically  pulsed  flow,  the  steady  state 
solution  at  a  given  streamwise  (?^)  station,  is  said  to  be 
reached  if 

f'(?,t  +  nT,n)  =  f'(?,t,n)  for  n  *  1.2,... 

where  T  is  the  period  of  oscillation. 

This  entails  velocity  profiles  to  agree  over  a  few  cycles 
and  thus  requires  unnecessary  and  uneconomic  computations 
extending  for  a  few  periods  in  order  to  ascertain  that  steady 
state  has  been  attained.  However,  if  all  the  three  para¬ 
meters,  f,  f'  and  f"  are  considered,  it  is  adequate  to 
regard  that  the  steady  state  solution  has  been  reached  if  at 
some  point  ns,  the  following  criteria  are  satisfied  - 

I  [f  U,tn  +  T  ,ng)  *  f(C,tn,ns)]/(f(^ ,tn,ns)  |  <  e4 

|[f’U,tn  +  T,ns)  -  f '  U,tn,ns)  ]/f '  (c ,tn,ns) !  <  e  n  =4?^2 

|[f”(?,tn  ♦  T,ns)  -  f"C5,tn,ns)]/f"(?,tn,ns)|  <  £(. 
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where  the  values  of  £^,  £,.  and  are  prescribed  and 
t2  -  t,  *  At  . 

The  sensitivity  of  the  steady-state  test  varies  with 
the  point  ns  to  which  it  is  applied.  In  practice,  the 
point  ris  should  be  so  chosen  such  that  reasonable  sensi¬ 
tivity  is  achieved  without  excessive  computations. 

Since  f’  and  f"  are  the  first  and  second  order  deriva¬ 
tives  of  f  ,  discrepancies  between  results  separated  by 
one  period  T  are  more  pronounced  in  f”  .  Hence,  in 
general  e4,  e5  and  must  be  so  chosen  that  e4  <  £,-  <  £&  . 
Otherwise,  if  e^  >  £^>  and  if  Eg  is  very  small,  the 
criteria  in  Eq.  (40)  might  never  be  satisfied  although  the 
steady  state  solution  has  long  been  attained  within  practical 
limits . 

As  a  further  check  of  the  attainment  of  the  steady-state 
condition,  the  following  condition  is  imposed  on  f’  at  the 
centre-line  - 

[[f'U,tn  +  T , 0 )  -  f’(C,tn,0)l/f'(C,tn,0)]  <  c7  (41) 
for  n  =  1,2 

where  the  value  of  Ey  is  prescribed. 

.  2 

It  has  been  chosen  that  *  Eg  =  Eg  =  Ey  =  10 
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3.7  Mean  Momentum  Flux,  Phase  Angle  and  Peak-to-Peak  Oscillation 
3.7.1  Mean  Momentum  Flux 

The  instantaneous  momentum  at  a  specified  streamwise 
station  is  given  by 


r  n„ 


Mi 


2hUjJ 


[f ' (s ,ti ,n) ] 1  dn 


(42) 


NL  is  evaluated  by  Simpson's  rule  for  unequally  spaced 
points  derived  in  Appendix  D  in  Lai  §  Simmons  (1978) . 

The  mean  momentum  flux  in  the  streamwise  direction  is 
defined  by 


M 


+  T 

dt 


(43) 


and  the  integral  is  approximated  by 


M 


N 

l 


(44) 


The  quantity  M  can  be  served  as  an  additional  check  on  the 
overall  results  by  testing  for  its  constancy  with  ?  . 

3.7,2.  Phase  Angle 

The  phase  angle  between  the  fundamental  component  of 
the  centre-line  velocity  at  any  downstream  station  and  that 
at  the  nozzle  can  readily  be  obtained  by  cross -correlating 
the  steady-state  instantaneous  centre-line  velocity  at  that 
station  with  a  reference  sine  and  cosine  signal  respectively. 
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Consider 

(45) 

(46) 


=  sin  cot 
S2  =  cos  cot 


From  Eq.  (39)  , 

f'c  =  f*(?,t,0)  -  fg  (?  ,t  ,0)  +  1  £a(C,  0)  sinUco  +  <f> , ) 

(47) 


Multiplying  Eq.  (45)  with  Eq.  (47)  and  taking  time  average 
yields 

3rj"T£  =  A  cos  cfrj  (48) 

where  A  is  a  constant. 

Multiplying  Eq.  (46)  with  Eq.  (47)  and  taking  time  average 
gives 

S,  =  A  sin  d>1  (49) 


The  percentage  peak  to  peak  variation  of  the  centre-line 
velocity  at  any  downstream  station  is  defined  by 


(£')  -  (f ' ) 

v  c^max  cJ 


x  1001 
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The  quantities  (f^max  and  *-f<Pmin  are  obtained  bY 
quardratic  interpolation  of  the  form 

f'c  =  aQ  t2  ♦  ax  t  +  a2  (52) 

through  three  points  (f^  (t^ ,  tj^) ,  (f^  (t2) ,  t2)  and 

(fc  V 

where  for  (f’)max,  f’^)  <  f'  Ct2) 

£'(t3)  <  f’  (t2) 

and  for  (f')nin  f’Ct,)  *  f’  (t2) 

f’(t3)  >  f’  (t2) 

Eq .  (52)  can  be  solved  for  a Qt  a^  and  a2  by  substituting 

f-(tl),  f^^)  and  ao»  ai  and  a2  known> 

^fc-*max  or  ^c^min  can  obtained  from  Eq.  (52)  with  t 
given  by  -a1/2a0  . 
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4.0  Results  and  Discussions 
The  structure  and  listing  of  the  computer  programs  for 
both  the  steady  and  unsteady  jet  are  described  in  Appendix 
A.  Calculations  were  performed  on  the  University  of 
Queensland  PDP105S  and  the  U.S.  Naval  Postgraduate  School 
IBM  360/67  computers.  The  computer  program  for  the  steady 
jet  occupies  a  core  memory  of  16K  words  and  that  of  the 
unsteady  jet  occupies  a  core  memory  of  67K  words.  It  must 
be  pointed  out  here  that  the  unsteady  jet  program  can 
actually  be  reduced  to  about  25K  words  of  core  memory  since 
at  any  one  time  instant,  only  calculations  involving  two 
time  levels  are  required  in  core  whereas  all  the  quantities 
at  other  time  levels  can  be  stored  on  disk  and  retrieved  when 
required.  However,  in  this  case,  the  computer  core  storage 
is  not  a  problem  with  the  computing  facility  available 
whereas  the  computing  time  is  important  and  as  such  a  trade¬ 
off  is  made  such  that  quantities  at  all  time  levels  are 
retained  in  core  in  order  to  save  computing  time  of  writing 
to  and  reading  from  disk.  Moreover,  although  a  constant  eddy 
viscosity  model  was  used,  the  programs  were  written  to  accept 
variable  eddy  viscosity.  A  very  straightforward  modification 
of  the  program  by  incorporating  a  routine  to  read  from  and 
write  to  a  disk  can  reduce  its  size  to  enable  it  to  be  run 
on  any  mini-computer  with  available  core  memory  of  about  25K 
words . 


4.1  Steady  Jet 

Owing  to  the  nature  of  the  initially  top-hat  velocity 
profile,  the  initial  velocity  gradient  is  very  steep  and  the 
velocity  profile  changes  vary  appreciably  over  a  very  short 
strearawise  distance.  Since  the  solution  of  the  finite 
difference  equations  is  obtained  through  linearization  by 
Newton's  method  as  described  in  section  3.3,  an  initial  gues 
at  any  streamwise  station  must  be  close  to  the  solution. 
Consequently,  a  fine  grid  has  also  to  be  used  initially  but 
as  the  computation  proceeds  and  as  the  velocity  profiles 
start  to  appear  in  similar  form,  a  coarser  grid  can  be  used 
by  dropping  every  other  point.  The  initial  grid  used  in 
the  transverse  (n)  direction  is  specified  as  follows: 


An  = 


CnJi 


0.005 

0 

< 

n 

< 

0.01 

0.002 

0.01 

< 

n 

< 

0.036 

0.0002 

0.036 

< 

n 

< 

0.039 

0.0001 

0.039 

< 

n 

< 

0.04 

0.00025 

0.04 

< 

n 

< 

0.045 

0.001 

0.045 

< 

n 

< 

0.06 

0.002 

0.06 

< 

n 

< 

0.02 

0.01 

n 

< 

0.12 

0.04 

(53) 


The  grid-sizes  in  the  streamwise  (?)  direction  are 
specified  as  follows: 
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0.0005 


0.005 

0.025 

0.1 

A?  = 

1 

2 

5 


0  <_  ?  _<  0.007 
0.007  <  5  <  0.032 
0.032  <  5  <  0.132 
0.132  <_  ?  <  1.032 
1.032  <  ?  <  20.032 
20.032  <  5  <  40-032 
40.032  <  5 


(54) 


The  constants  and  c2  in  equation  (12)  in  the 

constant  eddy  viscosity  model  are  varied  to  match  the  experi¬ 
mental  data.  The  sensitivity  of  the  solutions  to  various 
values  of  and  c2  was  tested.  No  significant  difference 

was  found  when  c^  was  varied  from  0.009  to  0.012  and  c2 
from  0.032  to  0.037.  The  final  values  of  *  0.009  and 
c2  ■  0.034  were  chosen. 

As  shown  in  Fig.  3,  the  non-dimensional  self -preserved 
velocity  profiles  of  various  workers  agree  very  well  with 
each  other  and  with  the  Goertler  solution  except  near  the 
jet  edge  where  conventional  hot-wire  measurements  are  dubious. 
However,  the  results  of  Heskestad  (1965)  follow  very  closely 
the  Goertler  solution  and  those  of  Robins  (1971)  which  have 
been  recommended  by  Rodi  (1975)  as  reliable  agree  very  well 
with  the  results  of  Heskestad  (1965).  The  computed  develop¬ 
ment  of  the  non-dimensional  velocity  profiles  of  the  steady 
jet  with  y  normalized  with  respect  to  the  jet  half-width 
y^  are  depicted  in  Fig.  4  and  at  streamwise  station  ?  * 
22.032,  the  mean  velocity  profile  differs  insignificantly  from 
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the  Goertler  solution,  indicating  that  self-preservation  in 
mean  velocity  profile  is  attained. 

The  variation  of  the  non-dimensional  centre-line  velocity 
with  streamwise  distance  is  plotted  in  Fig.  5.  The  results 
agree  well  with  the  experimental  data  of  Lai  and  Simmons 
(1979)  and  Zijnen  (1958).  The  variation  of  the  jet  half¬ 
width  with  streamwise  distance  is  compared  with  the  result 
of  Kotsovinos  (1976)  and  Lai  and  Simmons  (1979)  in  Fig.  6  and 
is  found  to  give  reasonably  good  agreement  for  streamwise 
distance  less  than  about  40  nozzle  widths.  The  discrepancy 
which  exists  at  larger  streamwise  distance  is  not  as  significant 
as  it  appears  because  the  results  of  Kotsovinos  (1976)  were 
obtained  by  fitting  a  third-order  polynomial  through  the 
scattered  data  in  the  literature  and  most  of  the  experimental 
data  were  not  available  for  streamwise  distance  larger  than  50 
nozzle  widths.  Furthermore,  uncertainity  in  the  hot-wire  data 
increases  as  the  streamwise  distance  increases  because  of  the 
decay  of  the  velocity  profile  and  possible  three-dimensional 
and  reversed  flow  effects.  The  rate  of  decay  of  the  centre- 
line  velocity  d  Uc/d(x/h)  of  0.165  and  the  spreading 
rate  d(yJj/h)/d(  x/h )  of  0.106  obtained  from  Figs.  5  and 
6  respectively  fall  within  the  range  of  values  reported  in  the 
literature.  As  pointed  out  by  Goldschmidt  and  Bradshaw  (1980), 
exact  self-preservation  in  mean  quantities  requires  the 
kinematic  virtual  origin,  obtained  from  the  centre-line 
velocity  decay  curve,  be  equal  to  the  geometric  virtual  origin 
obtained  from  the  jet  spreading  curve.  However,  most 
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experimental  data  do  not  confirm  this.  In  this  study,  the 
computed  kinematic  virtual  origin  agrees  with  the  geometric 
virtual  origin  and  is  found  to  be  0.8h. 

The  development  of  the  normalized  shear  stress  profiles 
is  shown  in  Fig.  7.  The  agreement  between  the  computed  shear 
stress  profile  at  ?  =  95. 032  and  the  experimental  results 
of  Gutraark  and  Wygnanski  (1976)  is  good  considering  that  the 
shear  stress  term  is  the  most  difficult  to  be  measured  with 
sufficient  accuracy  and  a  scatter  of  at  least  20%  exists  in  the 
available  experimental  data  in  the  literature.  The  maximum 
value  of  the  computed  non-dimensional  shear  stress  term 
u ' v* /U c  is  0.023  which  agrees  well  with  most  measured  values 
reported  in  Rodi  (1975). 

In  jet  computations,  numerical  errors  introduced  at  the 
jet  edge  normally  have  little  overall  influence  on  the  calcu¬ 
lations  as  noted  by  McGuirk  and  Rodi  (1979).  The  momentum 
integral  in  this  study  varies  by  less  than  1%  over  a 
streamwise  distance  of  100  nozzle  widths. 
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4.2  Unsteady  Jet 

The  constant  eddy-viscosity  model  which  was  found  to 
give  good  agreement  with  experimental  data  was  used  in  the 
computation  of  the  preiodically  pulsed  jet.  Results  were 
obtained  for  three  frequencies  of  pulsation,  u>  =  0.000871, 
0.00871  and  0.0871  which  for  a  jet  exit  Reynolds  number  of 
10  and  a  jet  width  h  of  5om  correspond  to  1 ,  10  and 
100  Hz  respectively.  Two  values  of  amplitude  of  pulsation 
were  studied,  namely  e  *  0.1  and  0.15.  The  results  were 
compared  with  the  experimental  data  of  Lai  and  Simmons  (1979) 
and  solutions  for  steady  jets.  The  grid-sizes  used  in  the 
streamwise  (?)  and  transverse  (n)  directions  are  the  same  as 
those  of  the  steady  jet  given  in  section  4.1. 

4.2.1  Sensitivity  of  the  Solution  to  the  Convergence  Criterion 

Results  show  that  solutions  are  very  sensitive  to  the  con¬ 
vergence  criterion  applied  to  points  very  near  to  the  edge  of 
the  jet.  Normally,  the  closer  the  point  ncon  is  to  the  edge 
of  the  jet,  the  more  difficult  it  will  be  for  the  solution  to 
converge.  It  was  found  by  Lai  and  Simmons  (1979)  that  the 
convergence  criterion  can  be  applied  to  a  point  in  the  jet 
which  varies  with  streamwise  distance  and  is  defined  by 


^con 


(Oi  '  3  for  ?  <  0.1 

6(n«)min/10'1  for  ?  *  0<1 


(55) 


This  allows  a  comparatively  sensitive  measure  of  the 
convergence  of  the  solution  and  yet  maintains  a  reasonably 


fast  convergence  without  affecting  the  overall  accuracy  of 
the  solution. 

A  fixed  point  arbitrarily  defined  by  ncon  =  r\m  -  6  may 
also  be  sufficient. 

4.2.2  Validity  of  the  Criterion  for  the  Attainment  of  Steady 
State  Solution 

The  criterion  for  the  attainment  of  steady  state  solu¬ 
tion  described  in  section  3.6  is  found  to  be  adequate.  The 
point  ns  to  which  the  steady  state  tests  are  applied  has 
been  chosen  to  coincide  with  ncon  specified  in  equation  (53). 

4.2.3  Sensitivity  of  the  Solution  to  Time  Step 

The  sensitivity  of  the  solution  to  time  step  has  been 

tested  for  various  frequencies.  Calculations  have  been  carried 

1 

out  with  the  period  of  pulsation  divided  into  12,  and  49  inter¬ 
vals  respectively .  Results  indicate  that  although  division 
of  a  period  into  12  time  intervals  seems  to  be  too  coarse,  they 
differ  insignificantly  from  each  other.  It  is  therefore 
sufficient  to  use  12  time  intervals  in  the  present  study. 

4.2.4  Relaxation  of  the  Criterion  for  Jet  Spreading 

As  pointed  out  in  section  3.5,  the  jet  spreads  rapidly 
initially  due  to  a  steep  velocity  gradient.  As  the  computa¬ 
tions  proceed  downstream,  the  velocity  profile  becomes  wider 
with  a  long  tail  at  the  jet  edge.  Because  of  the  very  low 
velocities  at  the  jet  edge,  small  numerical  errors  introduced 
there  might  cause  instabilities  and  the  criteria  set  out  in 
equation  (37)  for  the  definition  of  the  jet  edge  will  not  be 
satisfied  resulting  in  continuous  addition  of  unnecessary  grid 
points.  Futhermore,  the  number  n  of  grid  points  which  can 
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be  added  each  time  imposes  a  trade-off  between  the  accuracy 
of  the  solution  and  the  computing  time  required.  This  is 
because  less  computing  time  will  be  required  if  more  points 
are  added  but  since  the  quantities  at  the  new  points  are 
obtained  by  extrapolation  and  the  jet  edge  is  particularly 
sensitive  to  small  numerical  errors,  instability  or  inaccurate 
solution  may  result.  The  program  provides  choices  of  two 
different  number  of  added  grid  points  as  the  calculations  are 
marched  downstream  and  if  the  jet  spreads  beyond  a  certain  n 
value  which  can  be  obtained  from  the  Goertler  solution,  a  new 
definition  for  the  jet  edge  is  applied  to  terminate  the  jet 
at  the  point  where  the  velocity  starts  increasing  or  changes 
sign. 

4.2.5  Results 

Non-dimensional  mean  and  instantaneous  velocity  profiles 
at  various  time  delay  intervals  are  plotted  in  Fig.  8  for  the 
various  frequencies  and  amplitudes  of  pulsation  and  ?  »  40.032. 

The  mean  velocity  profiles  follow  closely  the  steady  jet 
curve.  The  instantaneous  profiles  also  collapse  into  the 
steady  jet  curve  except  in  Figs.  8(e)  and  (f)  where  a  slight  dis- 
crep«ncy  occurs  for  u>  =  0.00871  and  0.0871  and  e  *  0.15.  This 
agrees  with  the  trend  observed  in  the  experimental  data  of 
Lai  and  Simmons  (1979),  in  which  a  slight  departure  of  the 
non-dimensional  instantaneous  velocity  profiles  from  the 
steady  jet  curve  exists. 

The  mean  centre-line  velocity  Uc  obtained  by  time¬ 
averaging  the  steady  state  instantaneous  centre-line  velocity 
Uc  over  a  period  is  shown  in  Fig.  9  to  agree  well,  with  the  centre-line  velocity 


decay  curve  for  the  steady  jets  for  all  the  computed  cases. 

Fig.  10  shows  the  variation  of  the  jet  half-width  of  the 
mean  velocity  profile  with  streamwise  distance.  For  all 
cases,  the  mean  jet  half -width  of  the  unsteady  jet  collapses 
on  the  steady  jet  spreading  curve.  A  plot  of  Q/QE,  a 
measure  of  mean  entrainment,  versus  £  for  various  frequencies 
and  amplitudes  of  pulsation  indicates  that  the  mean  entrainment 
rate  does  not  differ  from  that  of  the  steady  jet,  which  is 
consistent  with  the  trend  associated  with  the  mean  centre¬ 
line  velocity  decay  and  jet  spreading.  In  all  the  computed 
cases,  the  mean  momentum  flux  M  is  conserved  to  within  1% 
over  a  streamwise  distance  of  100  nozzle  widths,  thus 
confirming  the  accuracy  of  the  solutions. 

With  reference  to  Fig.  12,  the  non-dimensional  mean 
shear  stress  profiles  of  the  unsteady  jet  at  c  =  40.032 
agrees  very  well  with  that  of  the  steady  jet.  The  non- 
dinensional  instantaneous  shear  stress  profiles  can  also  be 
collapsed  into  the  steady  jet  profile  except  for  to  =  0.0871 
in  Figs.  12(c)  and  (f)  where  significant  discrepancies  are 
noted.  This  is  consistent  with  the  behavior  of  the  instan¬ 
taneous  velocity  profiles  where  a  slight  disagreement  will  be 
amplified  in  the  velocity  gradient  term  which  is  shown  up  in 
the  shear  stress  profile. 

Fig.  13  shows  the  variation  of  the  steady-state 
instantaneous  centre-line  velocity  with  time  for  various 
frequencies  and  amplitudes  of  pulsation  at  three  different 
streamwise  distances.  The  reference  signal  is  given  by 
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IT  (1  +  esin  wt) .  It  can  be  noted  that  as  the  frequency  and 
amplitude  increases  and  for  large  streamwise  distance,  distor¬ 
tions  in  the  waveform  of  the  centre-line  velocity  exist,  an 
indication  of  the  importance  of  higher  harmonics.  Such  a 
trend  was  also  observed  by  Lai  and  Simmons  (1979)  in  their  experi¬ 
mentally  pulsed  jet.  The  variation  of  the  phase  angle  $ 
between  the  steady-state  fundamental  component  of  the  centre¬ 
line  velocity  at  any  streamwise  station  and  that  at  the  nozzle 
exit  with  ?  is  depicted  in  Fig.  14.  The  phase  angle  is  a  lag 
which  increases  w,ith  both  frequency  and  streamwise  distance. 
However,  the  amplitude  of  pulsation  does  not  have  any  significant 
effect  on  the  phase  lag  variation.  For  a  region  sufficiently 
close  to  the  nozzle,  the  flow  varies  with  time  in  a  quasi¬ 
steady  manner.  That  is,  the  solution  over  the  region  of 
interest  can  be  approximated  by  a  sequence  of  steady  jet 
solutions,  each  of  which  corresponds  to  the  instantaneous 
conditions  at  the  nozzle.  The  region  over  which  quasi¬ 
steady  approximations  can  be  applied  to  an  unsteady  jet  can 
readily  be  determined  by  using  the  phase  plot  in  Fig.  14. 

It  is  taken  arbitrarily  here  that  a  phase  lag  of  5  degrees 
defines  the  downstream  limit  to  quasi-steady  solutions. 

This  definition  enables  the  quasi-steady  approximations  to 
Uc  to  be  accurate  to  within  1%.  From  Fig.  14  it  is  apparent 
that  the  region  over  which  quasi-steady  approximations  can  be 
used  decreases  with  increasing  frequency.  The  quasi-steady 
region  for  the  computed  frequencies  and  for  e  =  0.1  or  0.15 
is  as  follows  - 
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«  =  0.000871 


0  <  ^  <  60 
“  =  0.00871  0  _<  C  <  10 

“  =  0.0871  0  <  z  <  1 

The  variation  of  the  percentage  peak-to-peak  oscillation, 
F,  of  the  steady-state  instantaneous  centre-line  velocity 
expressed  as  a  percentage  of  the  mean  value  with  streamwise 
distance  is  plotted  in  Fig.  15.  Although  for  co  =  0.0871, 

F  differs  quite  significantly  from  its  initial  value  of  20? 
or  30?,  it  varies  very  little  for  other  tested  frequencies  and 
amplitudes  of  pulsation.  In  the  experimental  data  of  Lai 
and  Simmons  (1979) ,  it  was  shown  that  for  frequencies  between 
1  and  10  H'z,  which  correspond  to  w  =  0.000871  and  0.00871  here, 
the  general  trend  indicates  that  F  increases  to  a  maximum  at 
10  nozzle  widths  and  drops  off  beyond  that  point  and  F 
increases  by  almost  40?  over  the  initial  value  even  for 
1Hz  case.  The  present  computations  using  the  constant  eddy 
viscosity  model  fail  to  predict  this  trend. 
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5.0  Concluding  Remarks 

A  method  which  employs  a  transformation  to  solve  the  thin 
shear  layer  equations  for  the  steady  and  unsteady,  two-dimen¬ 
sional  turbulent  free  let  issuing  into  stationary  air  has  been 

t 

presented..*  A  Prandtl  type  constant  eddy  viscosity  formula¬ 
tion  was  used  to  model  the  Reynolds  shear  stress  term.  The 
transformation  reduces  the  rate  of  spread  of  the  jet  and 
enables  the  use  of  variable  grid  sizes  over  predetermined 
regions  leading  to  accurate  predictions. 

In  the  steady  jet,  the  mean  velocity  profiles,  the 
mean  centre-line  velocity  decay  rate,  the  spreading  rate  of 
the  jet  and  the  normalized  shear  stress  profiles  all  show 
good  agreement  with  experimental  data. 

In  the  unsteady  jet,  the  mean  flow  characteristics,  such 
as  the  mean  velocity  and  shear  stress  profiles,  the  mean 
centre-line  velocity  decay,  the  mean  jet  spreading  and  entrain¬ 
ment  all  follow  closely  those  of  the  steady  jet  for  the 
range  of  tested  frequencies  and  amplitudes  of  pulsation.  This 
is  in  agreement  with  the  only  available  experimental  data  of 
this  type  which  are  available  in  Lai  $  Simmons  (1979) .  How¬ 
ever,  unsteady  effects  are  apparent  in  the  distortion  of  the 
waveform  for  the  variation  of  the  instantaneous  centre-line 
velocity  with  time  and  in  the  increase  in  the  phase  lag  with 
streamwise  distance.  For  high  frequencies  and  amplitude  of 
pulsation  or  at  large  streamwise  distance  from  the  nozzle, 
instantaneous  quantities  depart  from  quasi-steady  values. 

The  flow  region  over  which  quasi -steady  approximations  are 
applicable  has  been  established  for  various  frequencies. 
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f— 4 


Although  the  constant-eddy  viscosity  model  used  gives 
good  predictions  for  the  steady  jet  and  the  mean  quantities 
for  the  unsteady  jet,  it  fails  to  predict  accurately  the 
instantaneous  quantities  especially  the  precentage  peak  to 
peak  oscillation  of  the  steady-state  instantaneous  centre¬ 
line  velocity.  This  suggests  that  if  only  mean  flow 
quantities  are  required  and  if  the  frequencies  and  amplitudes 
of  pulsation  are  of  the  order  of  those  in  this  study,  the 
constant  eddy -viscosity  model  is  adequate.  However,  if 
instantaneous  quantities  are  critical  and  if  frequencies  and 
amplitudes  of  pulsation  are  an  order  of  magnitude  higher 
than  those  in  this  study,  a  more  refined  turbulence  model  has 
to  be  sought  and  tested. 
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Appendix  A  Structure  and  Listings  of  Computer  Program 
A.l  Sturcture  of  the  Steady  Jet  Program 

The  structure  of  the  input  data  is  as  follows: 

CARD  1  NRG,  (RG (I) ,  I  =  1,  NRG),  DEGR 

FORMAT  15,  6F10.6 

CARD  2  PROD  ,  CONST  1,  CONST  2 

FORMAT  3F10.6 

CARD  3  CO,  XO,  AA,  BB,  NPG ,  INC,  I FREQ,  NXI 

FORMAT  4F10.6,  414 

CARD  4  (NXTT(I),  1=2,  NXII) 

FORMAT  1014 

CARD  5  (DELX(I),  1=1, NXI),  X(l) 

FORMAT  8F10.6 

CARD  6  (DETAC(I),  I  =  1,8) 

FORMAT  8F10.6 

CARD  7  (VC (I),  I  =  1,7),  ETAE 

FORMAT  8F10.6 

CARD  8  El,  E2,  E3 

FORMAT  3F10.6 

CARD  9  OUTDSK’,  LN 

FORMAT  LI,  13 

CARD  10  ETAG 

FORMAT  F10.6 
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The  symbols  used  have  the  following  meaning: 


NRG 

No.  of  times  grid  points  have  to  be  rearranged 

RG(I) 

Upper  n  value  above  which  grid  points  will 

not  be  dropped 

DEGR 

Lower  n  value  below  which  grid  points  will 

not  be  dropped. 

PROD 

v  defined  in  section  2.1  (Jet  Exit  Reynolds 

Number) 

CONST  1, 
CONST  2 

Cj^  and  c2  defined  in  Eq.  (12) 

CO 

C0%/(ab) 

XO 

?o 

AA 

a  defined  in  Eq.  (lSa) 

BB 

b  defined  in  Eq.  (ISb) 

NPG 

nff,  No.  of  grid  points  to  be  added  if  the  jet 
spreads 

INC 

interval,  in  terms  of  the  number  of  transverse 
grid  points,  at  which  a  value  of  the  velocity 
profile  is  required 

I  FREQ 

interval,  in  terms  of  the  number  of  streamwise 
stations,  at  which  velocity  profile  is  required 

NX  I 

No.  of  different  grid  sizes  in  streamwise  direction 

NXTT(I) 

Streamwise  station  number  at  which  a  different  grid 
size  is  used. 

PE  LX  (I) 

Streamwise  grid  size 

X(l) 

Value  of  r  at  initial  streamwise  station  (zero  at 
nozzle  exit) 

DETAC(I) 

Transverse  grid  size  An 

VC  (I  ) 

Values  of  n  for  each  sub-region 

ETAE 

COi 

Convergence  limit  defined  in  Eq.  (36) 
Jet  edge  definition  defined  in  Eq.  (37) 


El 

E2,  E3 


49 


OUTDSK 

LN 

ETAG 


Assumes  logical  value  TRUE  if  solutions  are 
written  to  disk 

Device  number  of  writing  to  disk 

n  defined  in  Eq .  (28) 
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The  structure  of  the  input  data  is  as  follows: 


CARD  1 

NRG,  (RG(I),  I  =  1,  NRG),  DEGR 

FORMAT 

IS,  6F10.6 

CARD  2 

(NXG(I),  I  -  1,  NRG) 

FORMAT 

515 

CARD  3 

NPG ,  NPG2 ,  START,  LN,  NEW, 

NTRANS 

FORMAT 

715 

CARD  4 

CONST  1,  CONST  2 

FORMAT 

2F10.4 

CARD  5 

CO,  XO,  EPS,  OMG ,  NT,  INC 

FORMAT 

4F10.4,  215 

CARD  6 

IFREQ,  IFR,  I FA,  NXI 

FORMAT 

415 

CARD  7 

(NXTT(I) ,  1*2,  NXI I) 

FORMAT 

1015 

CARD  8 

(DELX(I) ,1*1,  NXI) ,  X(l) 

FORMAT 

8F10.4 

CARD  9 

CDETAC(I) ,  I  =  1,8) 

FORMAT 

8F10.4 

CARD  10 

CVC(I),  I  =  1.7),  GTAG 

FORMAT 

8F.10.4 

CARD  11 

El,  E2 ,  E3 ,  E4 ,  E5,  AA,  BB, 

PROD 

FORMAT 

8F10.4 

CARD  12 

ETAG 

FORMAT 

F10.4 

NPRINT 
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Symbols  which  appear  also  in  the  steady  jet  program 
have  the  same  meanings  here.  All  other  symbols  are  defined 
as  follows :  - 

NXG(I)  Streamwise  station  at  which  grid  has  to  be 

rearranged. 

NPG  2  Second  choice  of  number  of  transverse  grid 

points  to  be  added  if  jet  spreads  as 
discussed  in  section  4.2.4 

START  Streamwise  station  number  at  which  computations 

start 

NEW  gee  below 

NTRAN3  Streamwise  station  number  at  and  beyond  which 

constant  c^,  is  used. 

NPRINT  Streamwise  station  number  at  and  beyond  which 

velocity  profiles-are  required 

EPS  e,  amplitude  of  pulsation 

OMG  co,  angular  frequency  of  pulsation 

NT  Number  of  time  intervals  in  a  period  +  1 

IFR  interval,  in  terms  of  the  number  of  streamwise 

stations,  at  which  instantaneous  velocity  profiles 
are  required. 

IFA  interval,  in  terms  of  the  number  of  streamwise 

stations,  at  which  mean  quantities  are  required 

E4  and  £y  defined  in  Eq.  (40)  and  (41) 

E5  Eg  defined  in  Eq.  (40) 

Because  of  the  initially  steep  velocity  gradient,  the 
velocity  gradient  curve  behaves  erratically  at  the  region 
around  the  initial  jet  edge  which  constitutes  about  31  of  the 
total  jet -width.  To  eliminate  such  numerical  erratic 
behavior,  a  first-order  smoothing  function  can  be  applied  and 
is  supplied  through  the  Subroutine  SMOOTH.  The  use  of  the 


smoothing  function  does  not  affect  the  overall  results  and 
is  optional.  The  parameter  NEW  specifies  the  streamwise 
station  number  where  smoothing  function  starts  to  be  used. 

If  smoothing  function  is  not  required,  NEW  can  assume  the 
value  of  the  final  streamwise  station  number. 
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Appendix  B.  Figures 


Figure  1.  Cinfiguration  of  the  Plane  Jet 
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FIGURE  2.  Net  cube  for  the  difference  equations  for  1-q  (21). 
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Figure  3,  Non-Dimensional  Self-Preserved  Velocity  Profile  of  the  Steady  Jet 


Figure  4(b)  Non-Dimensional  Velocity  Distribution  of  the  Steady  Jet 
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Figure  6  Variation  of  the  Half-Width  of  the  Steady  Jet  with  Streamwise  Distance 
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Figure  7(b)  Nondimens ional  Shear  Stress  Profile  for  the  Steady  Jet 
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Figure  12(e)  Non-Dimensional  Shear  Stress  Profile 
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.Figure  12(f)  Non-Dimensional  Shear  Stress  Profile 
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Figure  13(a)  Variation  of  Instantaneous  Centre-Line  Velocity  with  Time  for 
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Variation  of  Instantaneous  Centre-Line  Velocity  with  Time  for  t  «  38.032 


Figure  13(d)  Variation  of  Instantaneous  Centre-Line  Veloc 
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Figure  13(e)  Variation  of  Instantaneous  Centre-Line  Velocity  with  Time  for  x,  »  38.032 
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Figure  13(f)  Variation  of  Instantaneous  Centre-Line  Velocity  with  time  For  c  -  60.032 


Figure  14(b)  Variation  of  Phase  Mg le  with 
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